/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// errCheck -- functions that check the consistency and/or compatibility of variables, input parameters etc.
// When consistency is violated, warning/error messages are issued.

/* program include files */
#include "dataStructs.h"

/* external data reference(s) */
extern Site *site;  
extern const VecR *siteCoords, *siteVelocity, *siteForce; 
extern const int applyThermostat, applyBarostat, nSites, stepPdb, stepNebr;
extern const VecI initHalfCellWat;
extern const VecR region;
extern const real *mass, rNebrShell;
extern const Prop boxVolume; 
extern real rCutMax;

void PrintMatrix( RMat );

char *errorMsg[] = {"", "bond snapped", "read checkpoint data",
		    "write checkpoint data", "copy buffer full", "empty event pool",
		    "message buffer full", "outside region", "read snap data",
		    "write snap data", "subdivision unfinished", "too many cells",
		    "too many copied sites", "too many layers", "too many levels",
		    "too many sites", "too many moved sites", 
		    "too many neighbors - increase input parameter 'nebrTabFac'",
		    "too many replicas",  "undetermined pair interaction type",
		    "#sites built != input parameter nSites"};

void ErrExit (int code) 
{
  printf ("Error: %s\n", errorMsg[code]);
  exit (0);
}

void CheckBaroThermoParameters( const real tauP, const real tauT, const real deltaT ) 
{
  real tauT_over_deltaT = tauT / deltaT; // (tauT / deltaT) ratio between T-coupling constant and integration timestep
  real tauP_over_deltaT = tauP / deltaT; // (tauP / deltaT) ratio between T-coupling constant and integration timestep
  
  if ( (applyBarostat) && (!applyThermostat) ) {   
    printf("\n**Error** in input parameters: applyBarostat requires applyThermostat.\n");
    printf("When controlling the pressure, the temperature must also be controlled (NPT conditions).\n");
    printf("However, it is possible to control the temperature but not the pressure (NVT conditions).\n");
    printf("Exiting...\n\n");
    exit(0);
  }
  
  if ( applyThermostat ) {
    if ( tauT_over_deltaT < 1 ) {
      printf("\n* Warning: input parameter tauT (%.f fs) is very small in relation to deltaT (%.f fs).\n", tauT*TIME_fs, deltaT*TIME_fs);
      printf("* The temperature-coupling might be too tight and lead to instabilities.\n");
      fflush( stdout );
    } else if ( tauT_over_deltaT > 999 ) {
      printf("\n* Warning: input parameter tauT (%.f fs) is very large in relation to deltaT (%.f fs).\n", tauT*TIME_fs, deltaT*TIME_fs);
      printf("* The temperature-coupling might be too weak; check the temperature (average and fluctuations).\n");
      fflush( stdout );
    }
  }

  if ( applyBarostat ) {
    if ( tauP_over_deltaT < 5 ) {
      printf("\n* Warning: input parameter tauP (%.f fs) is very small in relation to deltaT (%.f fs).\n", tauP*TIME_fs, deltaT*TIME_fs);
      printf("* The pressure-coupling might be too tight and lead to instabilities.\n");
      fflush( stdout );
    } else if ( tauP_over_deltaT > 999 ) {
      printf("\n* Warning: input parameter tauP (%.f fs) is very large in relation to deltaT (%.f fs).\n", tauP*TIME_fs, deltaT*TIME_fs);
      printf("* The pressure-coupling might be too weak; check the pressure (average and fluctuations).\n");
      fflush( stdout );
    }
  }
}

void CheckMinimumImageConvention( real rMax ) 
{
  if ( 0.5 * Min3( region.x, region.y, region.z ) < rMax + rNebrShell ) { 
    printf( "\n*** ERROR: The maximum cutoff size is not consistent with the minimum image convention (MIC).\n" );  
    printf( "Cutoff sizes must be smaller than half the shortest size of the simulation region:\n" );
    printf( "Since (rCutoffMax+rNebrShell=%.2f nm) is larger than ", rMax + rNebrShell );
    printf( "half the shortest edge (0.5*min(region.x,region.y,region.z)=%.2f nm), ", 0.5 * Min3( region.x, region.y, region.z ) );
    printf( "the MIC is violated.\nSolution(s): either increase the (shortest) size of the simulation region or reduce the cutoff sizes.\n");
    exit(0); 
  }
}

void CheckNLipids( const int nLipidsPerEdgeXY, const int nLipids ) 
{
  if ( Sqr( nLipidsPerEdgeXY ) * 2 != nLipids ) {
    printf("Error in InitLipidBilayer(): the requested number of lipids is not valid! Exiting...\n");
    printf("The total number of lipids nLipids should be such that sqrt( nLipids / 2 ) is an integer. Exiting...\n");
    exit(0);
  }
}

void CheckPressureTensor( const RMat ijPressure )
{
  int n;
  RMat nPressure; // pressure tensor
  RMat kinEnDensTensor; // kinetic energy density tensor 
  RMat globVirTensor; // global virial tensor 

  MZero( nPressure.u ); /* initialize */

  printf("\nPressure tensor from ij pairs:\n");
  PrintMatrix( ijPressure );

  printf("Pressure tensor from n sites:\n");
  DO_SITE {
    VOuter( kinEnDensTensor.u, siteVelocity[ n ], siteVelocity[ n ] );
    MScale( kinEnDensTensor.u, mass[ site[ n ].type ] );
    MMAdd( nPressure.u, kinEnDensTensor.u ); /* accumulate kinetic contribution to the pressure */

    VOuter( globVirTensor.u, siteCoords[ n ], siteForce[ n ] );
    MMAdd( nPressure.u, globVirTensor.u ); /* accumulate kinetic contribution to the pressure */
  }
  MScale( nPressure.u, 1. / boxVolume.val );
  PrintMatrix( nPressure );
}

void CheckStepParameters( const int nLipids ) 
{
  if ( 0 == stepPdb || stepPdb < 0 ) { 
    printf("Error: input parameter 'stepPdb' is not set properly - refer to manual. Exiting...\n"); 
    exit(0); 
  }
  if ( 0 == stepNebr || stepNebr < 0 ) { 
    printf("Error: input parameter 'stepNebr' is not set properly - refer to manual. Exiting...\n"); 
    exit(0); 
  }
}

void CheckStructTypeConsistency( int totNSites, int nSites, int nPointMasses, 
				 int nSymmRigBodies, int nWaters ) 
{
/*   printf( "pointMasses: %d\n", nPointMasses ); */
/*   printf( "#symmetricBodies: %d\n", nSymmRigBodies ); */
  if ( totNSites != nSites || totNSites != nPointMasses + nSymmRigBodies ) { 
    ErrExit( ERR_STRUCT_BUILT_INCONSISTENT_WITH_INPUT );
  }
}

/****************************************************************************************************************
 * CheckWatInputParameters -- checks that the bulk water configuration from "initHalfCellWat" is consistent with 
 *                            the actual number of water molecules 
 ****************************************************************************************************************/
void CheckWatInputParameters( const int nWaters ) 
{ 
  int halfHydratingWaters = VProd( initHalfCellWat ); // #waters to hydrate one lipid monolayer (either at the top or bottom of region)  
  if ( ( nWaters / 2 ) != halfHydratingWaters ) { // if the # of hydrating water is not consistent with the water input parameter
    printf("*** ERROR: the input parameter nWaters should be equal to 2*(initHalfCellWat.x * initHalfCellWat.y * initHalfCellWat.z)\n" );
    exit(0);
  }
}
